;; math/random.scm - (c) rohan drape, 2000-2005

;; Evaluates to a uniformly distributed random number in [0,1) with a
;; mean value of 0.5.

(define randomU base-random-uniform)

;; Evaluates to a linearly distributed random number in [0,1) with a
;; mean value of 0.2929.  The density function is given by 'f(x) = 2 *
;; (1 - x)'.

(define (random-linear) (min (randomU) (randomU)))

;; Evaluates to a linearly distributed random number in [0,1) with a
;; mean value of 0.6969.  The density function is given by 'f(x) = 2 *
;; (x - 1)'.

(define (random-inverse-linear) (max (randomU) (randomU)))

;; Evaluates to a triangularly distributed random number in [0,1) with
;; a mean value of 0.5.

(define (random-triangular) (* 0.5 (+ (randomU) (randomU))))

;; Evaluates to an exponentialy distributed random number with a mean
;; value of '0.69315 / l'.  There is no upper limit on the value
;; however there is only a one in one-thousand chance of generating a
;; number greater than '6.9078 / l'.  The density function is
;; given by: 'f(x) = l ^ (-l * x)'.

(define (random-exponential l)
  (let ((u (randomU)))
    (if (zero? u)
	(random-exponential l)
	(/ (- (log u)) l))))

;; Idiom over 'random-exponential' where all results are in [0,1).

(define (random-exponential-normal)
  (clipR= (random-exponential 6.9078) 1.0))

;; Evaluates to a bilinear exponentialy distributed random number with
;; a mean value of 0 and where half of the results lie between '+- 1 /
;; l'.  The denisty function is given by: 'f(x) = 0.5 * l * (e ^
;; (-l * |x|))'.

(define (random-bilinear-exponential l)
  (let ((u (* 2 (randomU))))
    (if (zero? u)
	(random-bilinear-exponential l)
	(* (if (> u 1) -1 1) (log (if (> u 1) (- 2 u) u))))))

;; Evaluates to a guassian distributed random number with a mean value
;; of `mu' and where 68.26% of values will occur within the interval
;; +-`sigma' and 99.75% within the interval +-(3 * `sigma').  The
;; density function is given by: 'f(x) = (1 / SQRT (2 * pi * sigma))
;; EXP - ((x - u) ^ 2) / (s * (sigma ^ 2))'.

(define (random-guassian sigma mu)
  (let* ((n 12)
	 (half-n (/ n 2))
	 (scale (/ 1 (sqrt (/ n 12)))))
    (let ((sum (apply + (map (lambda (ignored) (randomU)) (iota n)))))
      (+ (* sigma scale (- sum half-n)) mu))))

;; Idiom over 'random-guassian' where all results are in [0,1).

(define (random-guassian-normal)
  (clipU (random-guassian (/ 1.0 6.0) 0.5)))

;; Evaluates to a Cauchy-distributed random number with a mean value
;; of 0 where half of the results lie in the interval +-alpha, and
;; 99.9% fall within +-318.3alpha.  The density function is given by:
;; 'f(x) = alpha / (pi * (alpha ^ 2 + x ^ 2))'.

(define (random-cauchy alpha)
  (let ((u (randomU)))
    (if (= u 0.5)
	(random-cauchy alpha)
	(* alpha (tan (* u pi))))))

;; Evaluates to a beta-distributed random number in [0,1).

(define (random-beta a b)
  (let ((u1 (randomU)))
    (if (zero? u1)
	(random-beta a b)
	(let ((u2 (randomU)))
	  (if (zero? u2)
	      (random-beta a b)
	      (let* ((y1 (expt u1 (/ 1 a)))
		     (y2 (expt u2 (/ 1 b)))
		     (sum (+ y1 y2)))
		(if (> sum 1.0)
		    (random-beta a b)
		    (/ y1 sum))))))))

;; Evaluates to a weibull-distributed random number.

(define (random-weibull s t)
  (let ((u (randomU)))
    (if (or (zero? u)
	    (= u 1.0))
	(random-weibull s t)
	(let ((a (/ 1.0 (- 1.0 u))))
	  (* s (expt (log a) (/ 1.0 t)))))))

;; Evaluates to a weibull-distributed random number.

(define (random-poisson l)
  (let loop ((v (exp (- l)))
	     (u (randomU))
	     (n 0))
    (if (>= u v)
	(loop v (* u (randomU)) (+ n 1))
	n)))

;; The thread local parameter holding the procedure used to acquire a
;; random number in [0,1].  Initially this is set to 'randomU'.

(define randomN_ (make-parameter randomU))

;; The procedure to get a normal random value.

(define (randomN) ((randomN_)))

;; Evaluates to a random number in [0,`right') with a distribution
;; determined by `randomN'.

(define (random-zero-right right)
  (* (randomN) right))

;; Evaluates to a random number in [`left',`right') with a
;; distribution determined by `randomN'.

(define (random-left-right left right)
  (+ (random-zero-right (- right left)) left))

;; Evaluates to a number that lies within a deviation from `center'
;; calculated by multiplying `perturbation' by `center'.

(define (random-perturbation center perturbation)
  (let ((deviation (* center perturbation)))
    (random-left-right (- center deviation)
		       (+ center deviation))))

;; Evaluates to a random integer in [0,`right') with a distribution
;; determined by `randomN'.

(define (randomI-zero-right right)
  (floorE (random-zero-right right)))

;; Evaluates to a random integer in [`left',`right') with a distribution
;; determined by `randomN'.

(define (randomI-left-right left right)
  (+ left (randomI-zero-right (- right left))))

;; Evaluates to an integer that lies within a deviation from `center'
;; calculated by multiplying `perturbation' by `center'.

(define (randomI-perturbation randomN 
			      center perturbation)
  (let ((deviation (floorE (* center perturbation))))
    (randomI-left-right (- center deviation)
			(+ center deviation))))

;; Evaluates to '#t' or '#f'.

(define (randomB) (> 0.5 (randomN)))

;; Syntax for random branching, either `left' or `right' is evaluated.

(define-syntax random-branch
  (syntax-rules ()
    ((_ left right) (if (randomB) left right))))
